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The J 1 -J 2 Ising model in the square lattice in the presence of an external field is studied by two 
approaches: the Cluster Variation Method (CVM) and Monte Carlo simulations. The use of the 
CVM in the square approximation leads to the presence of a new equilibrium phase, not previously 
reported for this model: an Ising-nematic phase, which shows orientational order but not positional 
order, between the known stripes and disordered phases. Suitable order parameters are defined and 
the phase diagram of the model is obtained. Monte Carlo simulations are in qualitative agreement 
with the CVM results, giving support to the presence of the new Ising-nematic phase. Phase 
diagrams in the temperature-external field plane are obtained for selected values of the parameter 
k = J 2 /IJ 1 I which measures the relative strength of the competing interactions. From the CVM 
in the square approximation we obtain a line of second order transitions between the disordered 
and nematic phases, while the nematic-stripes phase transitions are found to be of first order. The 
Monte Carlo results suggest a line of second order nematic-disordered phase transitions in agreement 
with the CVM results. Regarding the stripes-nematic transitions, the present Monte Carlo results 
are not precise enough to reach definite conclusions about the nature of the transitions. 

PACS numbers: 05.50.+q,64.60.Cn,75.40.-s 


I. INTRODUCTION 

Competing interactions are common in many natu¬ 
ral and artificial systems, like the presence of conflict¬ 
ing ferromagnetic and antiferromagnetic interactions in 
frustrated magnetic systems as spin glasses [1] and ul- 
trathin magnetic films [2-4], as well as competition be¬ 
tween an attractive and a repulsive part in the interac¬ 
tion between atoms and molecules of complex fluids [5- 
11]. Competition between conflicting interactions is also 
relevant in mathematical optimization problems, when 
decisions have to be made where not all the constraints 
can be satisfied simultaneously [12]. Frustration, the in¬ 
ability of a system to statisfy all local constraints, is a 
unifying concept in many natural and artificial systems. 
Competing tendencies usually are responsible for com¬ 
plex behavior like slow relaxation to equilibrium, strong 
metastability and rough energy landscapes [13]. This 
makes the study of such systems both very interesting 
and challenging. One of the characteristic outcomes of 
the presence of competing interactions in a system is the 
emergence of heterogeneous structures as the equilibrium 
or low energy states, like stripes, bubbles, clusters, dis¬ 
ordered phases and anisotropic behavior. 

One of the simplest models with competing interac¬ 
tions is the J 1 -J 2 Ising model on the square lattice. This 
model is defined as a simple extension of the square lat¬ 
tice Ising model, in which besides the nearest-neighbor 
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(NN) ferromagnetic or attractive interaction Ji < 0, a 
next-nearest-neighbor (NNN) antiferromagnetic or repul¬ 
sive interaction J 2 > 0 is added. The ground state of the 
model depends on the relative intensity of the competing 
interactions n = J 2 /1 Ji |. For k < 1/2 it is ferromagnetic 
and for n > 1/2 it has a stripe structure of alternating 
up and down rows of spins. There is no exact solution for 
the thermodynamics of the model. At zero external field 
it has been studied by a variety of techniques like cluster 
mean field theory, transfer matrix and Monte Carlo sim¬ 
ulations [14-22], considering both ferromagnetic J\ < 0 
and antiferromagnetic J\ > 0 nearest-neighbor interac¬ 
tions. The nature of the thermal phase transition from 
the stripes to a disordered phase for k > 1/2 was con¬ 
troversial. In the most recent studies combining Monte 
Carlo simulations and a series of analytical techniques 
it has been established that the line of phase transi¬ 
tions in the temperature versus k plane is first order for 
1/2 < k < 0.67 and is continuous with Ashkin-Teller 
critical behavior for n > 0.67. The critical exponents 
change continuously in this regime between the 4-state 
Potts model behavior at k = 0.67 to standard Ising crit¬ 
icality for k —> 00 [20, 21]. 

In comparison with the zero field case, the model 
in an external field has received much less attention. 
Queiroz [23] and Yin et.al. [24] studied the case with both 
NN and NNN interactions of the antiferromagnetic type, 
using transfer-matrix methods in conjunction with finite- 
size scaling and conformal invariance in the first reference 
and large scale Monte Carlo simulations in the second 
one. For k = 1 the ground state is striped at small fields 
and a row-shifted phase appears for 4 < h < 8. The latter 
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state consists of alternating ferro and antiferromagneti- 
cally ordered rows (or columns), with the ferromagnetic 
ones parallel to the field. In both studies it was observed 
a reentrance in the boundary stripes-paramagnetic upon 
lowering the temperature at constant field. In Ref. 24 
it was argued that the reentrant behavior may be due to 
the appearance of row-shifted (2 x 2) clusters that help 
to sustain striped (2 x 1) order at low temperatures, even 
for moderately large magnetic fields. The nature of the 
phase transitions points to a weak universality scenario 
with exponents departing slightly from the standard Ising 
values. 

A natural question when dealing with stripe forming 
systems is the possibility of existence of an intermedi¬ 
ate nematic phase. A nematic phase in this context is 
characterized by the presence of orientational order but 
without translational or positional order [11, 25, 26]. In 
this sense there are broken symmetry phases but with 
an intermediate degree of symmetry, higher than the less 
symmetric stripe phases in which both orientational and 
positional orders are present. Nematic phases associ¬ 
ated with intermediate stripe-like order are present in 
many quasi-two-dimensional systems like ultratliin ferro¬ 
magnetic films[2, 25, 27] and electronic liquids in which 
they may be relevant to understand high temperature 
superconductivity [28-30]. The conditions under which 
a system can sustain a nematic phase of this kind are 
still not completely clear. Strong evidence for the ex¬ 
istence of such phases have been found in systems with 
isotropic competing interactions at different scales, e.g. 
when a short range ferromagnetic interaction competes 
with a long range antiferromagnetic one decaying with a 
power law of distance, like the dipolar interaction [31, 32]. 
In this particular case a nematic phase is present but 
only quasi-long-range nematic order develops in 2D. This 
quasi-nematic phase emerges by breaking a continuous 
0(2) symmetry, similar to the Kosterlitz-Thouless phase 
transition in the 2D XY model. In these kind of models, 
smectic phases are suppressed at finite temperatures due 
to the strong fluctuations of the order parameter. When 
an external magnetic field is applied, competing dipolar 
interactions give rise to new and interesting phases. At 
small fields stripe phases are still present, although the 
direction aligned with the field is favoured energetically 
and it gets wider as the field is risen until an instabil¬ 
ity leads to a first order phase transition to a bubble 
phase at a critical external field value. At still higher 
fields there is a second transition from the bubble to 
an homogeneously magnetized phase, a saturated para- 
magnet. Another salient feature of the field-temperature 
phase diagram of the system with dipolar interactions 
in a field is a strong reentrant behavior. This has been 
observed in beautiful experiments on ultrathin films of 
Fe/Cu(001) [33, 34] and also in mean field approxima¬ 
tions [35-37]. As for the presence of nematic phases in 
an external field, this is still an open question in dipolar 
or electronic systems. Compared to the behavior of the 
dipolar frustrated systems, the J 1 -J 2 model in the square 


lattice stands at the opposite side: it has a very simple 
stripe phase, with long range orientational and positional 
order, absent in models with long range isotropic inter¬ 
actions. The disordered-stripe phase transition in this 
model corresponds to the breaking of the symmetry 
of the square lattice to the Z 2 symmetry of the stripe 
phase. Then, the question we try to answer in this work 
is: is it possible for the J 1 -J 2 model in the square lattice 
to sustain an Ising-nematic phase? , i.e. a phase with 
orientational but not positional order ? To give even a 
partial answer to such question is always a considerable 
challenge. This is because the nematic order parame¬ 
ter in stripe forming systems amounts to compute cor¬ 
relation functions in different space directions searching 
for a breaking of isotropy characteristic of these phases 
[38, 39]. Then, the most common analytical approaches 
for a one-particle order parameter, namely mean field 
theory, fails at detecting nematic phases and one must 
go beyond naive MFT to an approximation which al¬ 
lows to compute anisotropic correlations. In this work 
we have studied the J 1 -J 2 model both with and without 
an external field by means of two approaches: the Cluster 
Variation Method (CVM) which is a cluster mean field 
theory, and Monte Carlo simulations. The CVM allows 
for a systematic improvement upon naive MFT by con¬ 
sidering clusters of particles of increasing size in an exact 
way in the partition function. To our knowledge, this 
is the first time the CVM is used with the specific aim 
of searching for anisotropic correlations, for which it is 
particularly well suited. In fact, the first step beyond the 
mean field approximation in a lattice is the two-site or 
pair approximation, also knwon as Bethe-Peierls approx¬ 
imation [40, 41]. This amounts to consider in a exact 
form all clusters with two sites. This approximation is 
known to predict correctly that d = 2 is the lower crit¬ 
ical dimension for the Ising ferromagnet. Nevertheless, 
because it does not distinguish any geometric or spatial 
features in the sum over pairs of sites, it is not able to cap¬ 
ture rotation symmetry-breaking, a distinctive feature of 
orientational phases. The next degree of approximation 
in the square lattice is the plaquette or square approx¬ 
imation, which considers exactly clusters of four sites, 
i.e. first and second neighbors. We will show that this 
is enough to capture the presence of nematic phases in 
models with competing first and second neighbor inter¬ 
actions. We did not find evidences of nematic phases in 
the J 1 -J 2 model at zero external field within the square 
approximation in the CVM, but a nematic phase appears 
when a field is present, both in the CVM approach and 
in Monte Carlo simulations. 


II. THE CLUSTER VARIATION METHOD 

We give here a very brief description of the Cluster 
Variation Method, focusing on quantities that will be 
useful in the calculation below. There is a large liter¬ 
ature on the technical aspects of the method and the 
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interested reader can refer to references [41-44] for more 
comprehensive discussions of the method and its poten¬ 
tial for computing variational approximations to the free 
energy of different systems. 

Consider the variational free energy 


for n = 1,2,..., IV—1, with the normalization conditions: 

Trp[ 1 \i) = 1, i = 1,2,..., N. (3) 

The cluster functions G and the cumulant functions 
c/( ra ) are defined as follows: 


F t = Tr(p t H) + k B T Tr (p t In p t ), (1) 


where Tr means a trace or a sum over all the relevant 
degrees of freedom of the Hamiltonian H , and pt is a 
trial density matrix which satisfies the normalization con¬ 
straint Tr p t = 1. A systematic way for obtaining varia¬ 
tional approximations to the free energy is to express it 
in terms of a cumulant expansion. Following the exposi¬ 
tion by Tanaka [41], consider the n-body reduced density 
matrix for an N body system: 


Pt " ) ( 1 7 2,..., n) = Tr n+1 p\ n+1) (1,2, ...,n,n+l) (2) 

I 


G {1) (i) =Tr 
G ( - 2 \i,j)=Tr 


p^(i) In p^ii) =g {1 \i), 
P?\i,j)knp ( ' 2 \i,j) 


= 9 {l) {i) + 9 (1 \j) + 9 (2 \hj), 
G {3 \i,j,k) = Tr \p[ 3) (i,j,k) Inp[ 3) (i,j,k) 


( 4 ) 

( 5 ) 


= 9 {1 \i) +9 {1 \j) +9 {1] (k) +9 {2 \i,j) 

+ 9 i2 \j,k) + g i2 \i,k) + g {3) (i,j,k), (6) 


and so on. The largest cluster fuction G^ (which corre¬ 
sponds to the N-body entropy term in the variational free 
energy) can be written as a sum over all the cumulant 
functions in the form: 


G W (1> ■ ■ ■ ,N) = Tr [pW(l, lnp t W (l, ...,N) 

= £s (1) (*)+£.9 (2) (*>j) + £ 9 i3 \^T k ) + --- + 9 [N) ( 


i<j 


i<j<k 


( 7 ) 


In this way the variational free energy can be written 
in terms of an expansion in cumulant functions: 


reduced density matrices pf 1 ^ can be written as [16]: 


F t =Tr 

Hp[ N \l,2,...,N) 

+ 

P^= 2- 

i + tffcCfc 





k 


k B T 


xy i)( 4) 


*<?■ 


( 2 ) 


(bj) + 


£ 9 

i<j<k 


(3) 


H-b 5 (JV) (1> ■ • ■! N) 


( 8 ) 


(9) 


where the sum runs over all subclusters with k sites 
within cluster n , Ofc = and the k-point cor¬ 


relation functions are defined by = Trokp\ 
variational parameters fk must satisfy: 


(fc) 


The 


This form of the variational free energy allows to obtain 
systematic approximations to the true free energy of a 
model system by considering a maximal size of cluster 
to be summed exactly in the partition function. This is 
called the parent cluster. This amounts to truncate the 
cumulant expansion to a given degree and optimizing the 
resultant expression with respect to the reduced density 
matrices which must satisfy the reducibility and normal¬ 
ization conditions (2) and (3) respectively. The reduced 
density matrices represent all the subclusters contained 
within the parent cluster. In some applications a conve¬ 
nient way of implementing the variational approximation 
is to parametrize the density matrices in terms of corre¬ 
lation functions and consider these as variational param¬ 
eters instead. 

In the case of a system with Ising spins {Si = ±1}, the 


m 

d( k 


= 0. 


( 10 ) 


A hierarchy of approximations to the free energy can be 
constructed in this way. The simplest one corresponds 
to the 1-point approximation for the density matrices, 
the usual mean field approximation. The 2-point ap¬ 
proximation is usually called Bethe-Peierls approxima¬ 
tion [40, 41]. As discussed in the Introduction, the 
pair approximation is not able to detect orientation- 
dependent features of the equilibrium phases. In this 
work, we implemented the 4-point approximation in the 
square lattice, which is able to capture the emergence of 
anisotropic nearest-neighbor correlations or spontaneous 
rotational symmetry breaking. 
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III. Ji-J 2 ISING MODEL IN THE 4-POINT 
(SQUARE) APPROXIMATION 

The J 1 -J 2 Ising model on the square lattice is defined 
by the Hamiltonian: 

n = j l Y J s x s y + j 2 J2 s * s y ( n ) 

(xy) {(xy)) x 

where {S/, = ±1 ,2 = 1... TV} are N Ising spin variables 
and h is an external field, (xy) denotes a sum over pairs 
of nearest-neighbors and ((xy)) a sum over pairs of next- 
nearest-neighbors. In this work we consider J\ < 0 and 
J 2 > 0 representing ferromagnetic NN and antiferromag¬ 
netic NNN interactions respectively. The competition 
ratio is defined by re = > 0 . 

At zero external field the ground state of the model is 
ferro (if J\ < 0) or antiferromagnetically ordered (if J\ > 
0 ) for re < 1/2 and striped or superantiferromagnetic if 
re > 1/2. In the stripe phase the system can adopt one 
of four possible configurations as shown in Figure 1. 

At zero external field the model has been extensively 
studied [14-22], considering both ferromagnetic Ji < 0 
and antiferromagnetic Ji > 0 NN interactions and anti¬ 
ferromagnetic J 2 > 0 NNN interactions. The nature of 
the thermal phase transition from the stripes to a disor¬ 
dered phase for re > 1/2 was controversial. In the most 
recent studies combining Monte Carlo simulations and a 


series of analytical techniques it has been established that 
the line of phase transitions in the temperature versus re 
plane is first order for 1/2 < re < 0.67 and is continuous 
with Ashkin-Teller critical behavior for re > 0.67. The 
critical exponents change continuously in this regime be¬ 
tween the 4-state Potts model behavior at re = 0.67 to 
standard Ising criticality for re —> 00 [20, 21]. In the 


FIG. 1. Sketch of the ground state configurations for the Ji- 
J 2 model without external field. Empty circles: Si = — 1; 
filled circles: Si = +1. 

square approximation, the CVM free energy correspond¬ 
ing to the Hamiltonian (11) is given by [16]: 




F = Ji^Tr (S x S y p( xy )) + J 2 E Tr { S xSyP((xy))) 

{xy) {(xy)) 


h^Tr ( S x p x ) + k B T 

X 


Ey (1) fa) 


+ E 5 (2) fa’ y) + E 9 (2 \x,y) + E 9 (3 \x,y,z) + E 9 {4 \x,y,z,w) 

(xy) ((xy >> [xyz] 


( 12 ) 


In the last equation, the sums over x, (xy), ((xy)), [xyz], 
denote sums over all sites, NN pairs, NNN pairs, 
clusters of three sites and squares respectively. 

Expressing the cunrulant functions g in terms of the 
cluster functions G the variational free energy reads: 


F = Ji^Tr (S x S y p {xy) ) + J 2 E Tr { S x S yP((xy))) 
{xy) {{xy)) 

- h^Tr(S xPx ) + k B T EG (1) W 

X X 


- E ° {2) fa’ y) + E g(4) fa’ 

{xy) 


(13) 


Following the general approach outlined in Section II it 
is convenient to write the variational free energy in terms 


of correlation functions given by: 

m x = Tr (S x p x ) 

Ixy — TT \$x$yP(xy)) 

Cxy — Tv x ^y P {{xy ))) 

kyxw — TV (^Sy S X S W p\y XW j j 
dxyzw — Tv ^S x SyS z S w px 5 (14) 

which are related to the reduced density matrices by: 

Px = 2 (1 4 “ 'W'xSx) 

P{xy) — (1 “b m x S x + 'TTlySy ^xy^x^y) 

P{{xy)) — ^ (1 “b TFlx^x "b 'WlySy H"~ C-xy^x^y) 
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= — (1 + m x S x + rriySy + m z S z + m w S w 
T Ixw^x^w 4~ Iwz^w^z “f~ Izy^z^y T Ixy^x^y 

F Cxz^x^z H - Cyw^y^w T kyxw^y^x^w 

“t" k xwz S x S w S z T k WZ yS w S z S y T k zyx S z S y S x 


F dxyzw^x^y^z^w)' 


(15) 


Substituting these definitions onto (13) we get the vari¬ 
ational free energy of the J 1 -J 2 model in the CVM square 
approximation[16]: 


f — Ji ^ ] i X y f J2 y ) Cxy h y ) 

{xy) {{xy)) x 


TO, 


k B T 


J2 Tr (Pxlogp x )-J2 Tr (p( xy )l°9P( xy)) 

x (xy) 


F 


Yz Tr (p^iogpiui 


(16) 


After computing the traces one is left with an expres¬ 
sion for the variational free energy in terms of a set of 
correlation functions representative of the approximation 
considered. The form of equation (16) makes clear that 
up to the pair approximation the two directions in the 
square lattice enter in a completely symmetric way, the 
different pairs of sites are decoupled. It is in the last term, 
when square plaquettes are considered, that the coupling 
between different directions in space can lead to novel be¬ 
havior. The minimization of the free energy is in general 
a difficult task. The state (stationarity) equations are 
given by (10). It is possible to compute the derivatives 
and try to solve the set of coupled nonlinear equations of 
state. Instead of that, we preferred to minimize the vari¬ 
ational free energy numerically for given sets of external 
parameters. k B = 1 was set for all calculations. 


IV. CVM RESULTS IN AN EXTERNAL FIELD 

For k > | and small magnetic fields the ground state 
is striped (2 x 1). When Ji is ferromagnetic and J 2 anti¬ 
ferromagnetic, the stripe order is eventually destroyed by 
the presence of an external field, and all the spins become 
aligned with the field at h c = ±2(Ji + 2 J 2 ). In the case 
both interactions are antiferromagnetic, for —4 J 2 < h < 
4 J 2 the ground state is still (2 x 1), while in the interval 
4J 2 < h < 4Ji + 4J 2 and —(4Ji + 4J 2 ) < h < —4J 2 
it becomes row shifted (2 x 2) [23, 24]. The latter state 
consists of alternating ferro and antiferromagnetically or¬ 
dered rows (or columns), with the ferromagnetic ones 
parallel to the field. For higher fields the equilibrium 
state of the system corresponds to a saturated paramag- 
net. 

Here we consider the system at finite h with ferromag¬ 
netic NN (Ji < 0) and antiferromagnetic NNN (J 2 > 0) 
interactions for two different competition ratios n = 1 


and k = 0.6 for which the ground state is striped. With 
the aim of searching for purely orientational nematic-like 
phases, i.e. phases without positional order, we mini¬ 
mized the CVM free energy of Eq. (16) for the param¬ 
eters (related to the elementary square defined in (12)): 
17l x — TTlw — Ixw — lyzi lxy> lwz> C, k yX w = k Z yxi 
k XW z = k wzy and d. This choice implies possible orien¬ 
tational order along the xy or vertical direction. Note 
that local magnetizations on horizontal NN sites are al¬ 
lowed to be different in sign and also in absolute value. 
Correspondingly, the NN correlation functions l rs may 
be different not only between the horizontal and ver¬ 
tical directions but also between the two vertical ones. 
With these choices the values of NNN correlations c and 
square correlations d are unique. The values of m r: l rs , 
c rs , k rs t and d with r, s,t = x,y,z,w that minimize the 
free energy (16) were calculated numerically for differ¬ 
ent reduced temperatures (T/| Ji|) and reduced external 
fields {h/\J\\) using the routine NMinimize of the soft¬ 
ware Mathematica [45]. In order to distinguish between 
translational (positional) order and orientational order 
we defined suitable order parameters; two positional or¬ 
der parameters: 


and 


M f 


2m x F 2 m w 
4 


(17) 


M s 


2 m x — 2 m w 
4 


(18) 


describing ferromagnetic and stripe orders respectively. 
Note that, in case m x and m w are both finite but with 
different absolute values, a mixed phase with stripe order 
on a ferromagnetic background can be possible. In this 
cases, we have classified these phases as stripe ones. The 
orientational order parameter is defined as: 


Q — ^ (.Ixy T Iwz 2 / xw)- ( 19 ) 

In this case orientational order is finite whenever hori¬ 
zontal NN correlations are different from vertical ones, 
i.e. when there is a breaking of Z 4 symmetry in the NN 
correlation functions. With these definitions all order 
parameters take values in the range (—1,1). The phase 
diagram in the h/ \ J\ \ — T/ \ J\ \ plane for k = 1 is shown in 
figure 2. For this value of k the ground state is striped. 
A first difference in this phase diagram with respect to 
those in references 23 and 24 is the absence of the row- 
shifted phase at large fields. This is due to the ferro¬ 
magnetic character of the NN interactions in the present 
work. The absence of the (2 x 2) phase is probably re¬ 
lated with the absence of reentrance of the stripe phase 
in this case, at variance with the results reported for the 
model with antiferromagnetic NN interactions [23, 24], 
as pointed out above. At h/\ J\ \ = 0, the 4-point approx¬ 
imation predicts a first order transition line for k < 1 
and a second order transition line for k > 1 [14 16]. For 
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FIG. 2. (Color online) Reduced external field versus reduced temperature phase diagram for k — 1. Full (dotted) lines 
correspond to continuous (discontinuous) transitions. 




h/IJjl h/IJjl h/IJjl 

FIG. 3. (Color online) Correlation functions versus reduced magnetic field for k = 1 and T/\Ji\ = 0.8. Panel a: positional 
(M s ) and orientational order parameters ( Q ). Panel b: local magnetizations. Panel c: nearest-neighbor correlations. Panel d: 
next-nearest-neighbor correlations. Panel e: three site correlations. Panel f: square correlations. 


k = 1 and finite h/\Ji\ we found a line of first order tran¬ 
sitions marked by the discotinuity of the order param¬ 
eter. The discontinuity is observed at both transitions, 
stripes-saturated paramagnetic and stripes-nematic (see 
Figure 3, first panel). Above the point where the posi¬ 
tional order parameter goes to zero the system still finds 
itself in a phase with a finite value of the orientational 
order parameter Q , as seen in Figs. 3 and 5. This is the 
signature of a nematic-like phase in which the magne¬ 


tization is homogeneous, unlike in the stripe phase, but 
correlations show an anisotropic character, reminiscent 
of the more ordered stripe phase. For the case k = 1 the 
nematic phase is observed in the h/ \ J\ | — Tj \ J\ | plane 
in the reduced temperature range 0.1 < T/\ J\\ < 1.4, as 
seen in Fig. 2. The nematic phase terminates in a line 
of second order phase transitions (green line in Fig. 2) 
where the system enters a paramagnetic phase with finite 
magnetization values due to the external field (saturated 
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paramagnet). The observation of the nematic phase in 
the J 1 -J 2 model in an external field is the main result of 
this work. In the context of the Cluster Variation Method 
it is clear that the 4-point approximation is the minimal 
one which is able to capture a nematic-like phase of this 
kind, i.e. a phase with broken orientational symmetry in 
the correlation functions. Nevertheless, this possibility 
was not exploited in previous work within the CVM. 

The behavior of correlation functions for k = 1 and 
T/\Ji\ = 0.8 as functions of the external field is shown 
in Figure 3. In the first panel of Figure 3 it is seen that 
the positional and orientational order parameters coin¬ 
cide in the stripe phase, as it should, but Ms goes to 
zero before Q , signalling a stripe-nematic phase transi¬ 
tion at h/\Ji\ ~ 1.96. The second panel on the upper row 
shows that the local magnetizations m x and m w tend to 
be equal but with opposite signs at very small fields, in 
agreement with the stripe character of the ground state 
for small h/\J±\ values. Nevertheless, they gradually 
evolve in an asymmetric way until at the stripe-nematic 
transition point their values merge in a single one, mean¬ 
ing the onset of an homogeneous phase with regard to 
local magnetizations. The anisotropic character of the 
nematic phase is evidenced in the third panel on the up¬ 
per row of Figures 3 and 5. It is seen that, within the 
stripe phase, correlations along the vertical direction ( l xy 
and l wz ) are slightly different reflecting the slightly differ¬ 
ent values of the local magnetizations at those sites. A 
change in behavior between those correlations and the 
horizontal ones l xw is observed at the stripes-nematic 
transition point. Note that if the stripe phase should 
have ended in a disordered rotationally symmetric state, 
then correlations in different directions should merge at 
this point. This does not happen until a larger field value 
where the three different NN correlations considered here 
merge to a single value at h/\J\\ ~ 2.08. The NNN cor¬ 
relations c display a discontinous derivative at the stripe- 
nematic transition. 

In Figures 4 and 5 we show the phase diagram and 
correlation functions for « = 0.6 and T/\Ji\ = 0.65 . 
The behavior is qualitatively the same as the case with 
k = 1. We show them in order to compare with Monte 
Carlo simulation results to be discussed in the next sec¬ 
tion for k = 0.6. Although the results from the CVM 
and Monte Carlo simulations are qualitatively similar, 
important quantitative differences arise which should be 
addressed with higher order approximations in the CVM 
approach. 


V. MONTE CARLO SIMULATIONS 

We have performed Monte Carlo simulations of the 
J1-J2 Ising model on the square lattice to test the qual¬ 
itative consistency of the CVM results. The frustration 
present in the model turns computer simulations very 
demanding. The simulation procedure combines stan¬ 
dard one-site moves [46] with a cluster method adapted 


from the simulation of patchy lattice models [47] follow¬ 
ing the ideas of the Wolff algorithm [48] in the presence 
of external fields [46]. In the cluster method, one of the 
spins of the system, and one of the main directions of 
the lattice are chosen at random. The chosen spin is 
the starting point (root) of the cluster, then one starts 
growing the cluster by adding spins, j, which are NN 
of the cluster in the chosen direction with probability 
b = 1 — exp(—2| Jil/ksT) provided that Sj = So, with 
So being the state of the spins in the cluster. Once the 
construction of the cluster is finished acceptance crite¬ 
ria are applied by taking into account the interaction of 
the cluster with the remaining spins of the lattice and 
the external field[46, 47]. The introduction of the clus¬ 
ter technique improves appreciably the numerical perfor¬ 
mance of the simulations, specially at low temperature. 
In order to enhance further the simulation peformance 
we have made use of parallel tempering (or replica ex¬ 
change) Monte Carlo sampling[49, 50]. This was carried 
out as follows: for a given fixed value of the external 
field (or temperature) we located via preliminary sim¬ 
ulations the approximate value(s) of the temperarature 
(or field) where the transition(s) take(s) place; then we 
chose an appropriate set of values of the temperature 
(field) around such preliminary estimates to carry out 
the replica-exchange Monte Carlo simulations. The ini¬ 
tial configurations were built by choosing the value of 
each spin at random. In order to guarantee the relia¬ 
bility of the final results we run shorter complementary 
simulations, using ground state configurations as starting 
point, to check that the simulations runs were properly 
equilibrated. 

In the same spirit as with the CVM approach, our 
main interest was to search for possible nematic phases, 
i.e. phases with orientational order but lacking trans¬ 
lational order. In order to distinguish the presence of 
such phases we defined suitable order parameters analo¬ 
gous to those defined previously in the CVM approach. 
The translational order parameter (TOP) was defined as 
the one used in Ref. 21, in which the configurations of 
the system are basically compared with the ground state 
configurations at zero field. We analyzed the existence of 
periodicity in the lattice directions a = x, y by comput¬ 
ing the quantities: 

1 N 

O t (a) = t 2 x mod (a*, 2) — l]Sn (20) 

i— 1 

where L is the linear size of the square lattice and 
oti = Xi, yi are the coordinates of site i. The global trans¬ 
lational order parameter Ot is then defined through: 

0 2 T = 0 2 t {x) + 0 2 t {y). (21) 

Notice that Ot = 1 for the ground state structures shown 
in Fig. 1. This order parameter is equivalent to the defi¬ 
nition used in (18) which was suitable for the elementary 
square of the CVM. The translational order can be tested 



FIG. 4. (Color online) Reduced external field versus reduced temperature phase diagram for k = 0.6. Full (dotted) lines 
correspond to continuous (discontinuous) transitions. 
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FIG. 5. (Color online) Correlation functions versus reduced magnetic field for k = 0.6 and T/\Ji\ = 0.65. Panel a: positional 
(M 3 ) and orientational order parameters ( Q ). Panel b: local magnetizations. Panel c: nearest-neighbor correlations. Panel d: 
next-nearest-neighbor correlations. Panel e: three site correlations. Panel f: square correlations. 


by carrying out a finite-size scaling analysis. If no trans¬ 
lational order exists one expects that the average of this 
order parameter (Ot) will approach zero in the thermo¬ 
dynamic limit, whereas for the ordered case (Ot) will 
be finite as L -> oo. Moreover, if we analyse the trans¬ 
lational order through an isotherm (at varying external 
field), or as a function of the temperature (at constant 
external field), we could expect that the possible order- 
disorder transitions involving translational ordering will 


appear as abrupt changes in (Ot), specially for large sys¬ 
tems. 

According to the ground states at zero field (Fig. 1), 
at low temperatures the system shows the tendency to 
form long sequences of spins in the same state along the 
main directions of the lattice x and/or y. The length 
of these sequences is expected to grow on decreasing the 
temperature and at some point the competition between 
sequences in the two directions might lead to a phase 
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transition, in such a way that one of the directions will 
be preferred in the ordered phase. In this case the system 
breaks the fourfold (Z 4 ) symmetry of the square lattice 
reducing it to twofold (Z 2 ) symmetry. This transition 
from a disordered Z 4 to a Z 2 ordered phase is an orien¬ 
tational phase transition and need not be accompanied 
with the growth of translational or positional order. So 
we must distinguish positional and orientational order 
parameters, as discussed in relation to the CVM results. 

We can define orientational order parameters (OOP) 
along directions x and y as: 

L 2 

°o{a) = ' 52 S(r i )S(r i + d); a = x,y\ ( 22 ) 

i=l 

We will find that for the GS configurations one of the 
components of O a will have the value +1 (that corre¬ 
sponding to the direction in which the sites are at the 
same state), whereas the other will take the value —1. A 
global order parameter is then defined as: 

( - > ° — \ \°o{x) ~ O 0 (y)\ ■ (23) 

This definition is equivalent to Eq. (19) which was suit¬ 
able in the context of the square approximation of the 
CVM. In the thermodynamic limit, finite values of the 
statistical average of its absoute value < \0 o \ > will indi¬ 
cate that there is a preferential direction for the sequences 
of equal spins. 


A. Results 

We have considered two values of k, namely k = 0.6, 
and n = 1.0. For n = 0.6 we have analysed five cases. 
In the first one, ft = 0 was fixed and we looked at the 
variation of different properties with the temperature and 
the system size. In addition, we took four values of the 
reduced temperature T/\Ji\ = 0.80, 0.6667, 0.50, and 
0.40, and looked at the variation of the properties as a 
function of the external field, ft, and the system size at 
constant temperature. 

Some results for k = 0.6, ft = 0 are shown in Fig. 6 . 
In this case the system exhibits a disordered phase at 
high temperatures, where both order parameters tend to 
zero as the system size increases. At reduced tempera¬ 
ture T c /\ Ji| ~ 0.972 a phase transition occurs. A direct 
inspection of the results (See Fig. 6 ) for different proper¬ 
ties as a function of T/ \ J\ | for different system sizes leads 
to the following conclusions: (1) For T > T c both order 
parameters seem to vanish as L grows larger; whereas 
for T < T c both order parameters converge to a finite 
value greater than zero for the larger system sizes; ( 2 ) 
The jump in both order parameters seems to occur at 
the same value of the temperature (T c ); (3) At T c there 
is also a jump in the energy per site: H/L 2 , (Fig. 6 .c); 
(4) The corresponding heat capacity at constant field, 


Ch = (d['H/L 2 ]/dT)h exhibits a clear single peak for sys¬ 
tems with L > 32, with a value for the maximum that 
seems to diverge as L —> 00 (Fig 6 .d). The scaling of the 
maximum of Ch for the system sizes considered does not 
allow to establish the type of the transition. Two possible 
scenarios for the transition are, a 4-state Potts criticality 
(continuous transitions), with c/, oc I 1 , or a very weak 
discontinuous transition as stated by Jin et al. [21]. 

Next, we consider the case k = 0.6, T/|Ji| = 0.8 and 
finite ft. The main difference with the case with ft = 0 is 
that the averaraged magnetization of the system, defined 

as (to) = L~ 2 (JA =1 Si), does not vanish in the thermo¬ 
dynamic limit. In Figure 7 we show some of the results 
for this system. From the results of the order parameters 
it can be deduced that at reduced field ft/|Ji| — —0.305 
the system exhibits an order-disorder transition. The re¬ 
sults suggest that the orientational and translational or¬ 
dering occurs simultaneously. At the same value of ft it is 
observed a jump in the magnetization, whose derivative 
with respect to the external field at constant tempera¬ 
ture, xt, seems to diverge as i A 00 . The qualitative 
behavior of the energy per site with respect to ft (not 
shown), and its derivative Ch as a function of ft is similar 
to the behavior for the case ft = 0 when both functions 
are plotted as functions of T. Therefore, the features 
of the transition at T/\Ji\ = 0.80 are similar to those 
found at zero field. In both cases the orientational and 
translational orderings seem to occur cooperatively. This 
conclusion seems to be consistent with the presence of an 
unique peak in the susceptibility \t, as shown in Fig. 7, 
panel d. 

On decreasing the temperature the qualitative features 
of the transitions change. We have carried out simula¬ 
tions at three additional values of reduced temperature: 
T/|Ji| = 0.6667, 0.50, and 0.40 (for k = 0.6). Some 
qualitative differences arise between the phase transitions 
in these three low-temperature cases and the preceding 
two cases. In Fig. 8 we show the results for the case 
T/|Ji| = 0.6667. First, the aparent common transition 
for TOP and OOP seems to split into two separated tran¬ 
sitions, i.e. for a given temperature the jumps of TOP 
and OOP occur at different values of the external field ft. 
This scenario of two successive order-disorder transitions 
instead of just one transition is consistent with the incip¬ 
ient splitting of the peak of the susceptibility \t for the 
larger systems considered. This behavior is in qualita¬ 
tive agreement with the CVM results for corresponding 
parameter values observed in Fig. 5. 

We have also explored the phase behavior for k = 1.0. 
Three cases were considered (1) Constant field ft = 0; (2) 
Constant temperature T/\Ji\ = 1.0, and (3) Constant 
temperature T/|Ji| = 0.50. The phase behavior is sim¬ 
ilar to that found for n = 0.60. At zero held, a single 
transition is found at T c /|Ji| ~ 2.08. For T/|Ji| = 1.0, 
the intermediate nematic phase (with only orientational 
order) does not appear, with the ordering transition oc¬ 
curring at |ft|/| Ji| ~ 1.88. However for T/\ Ji\ = 0.5 the 
Ising nematic phase seems to be stable for a very narrow 
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FIG. 6 . (Color online) Results for the J 1 -J 2 model with k = 0.60 and h = 0. Different system sizes (as indicated in the legends) 
are considered. The results are shown as a function of the reduced temperature. Panel a: Orientational order parameter 
(OOP); Panel b: Translational order parameter (TOP); Panel c: mean energy per site <T-L> /L 2 \ Panel d: Heat capacity per 
site at constant held Ch = (d[H/L 2 ]/dT)h- 


range of h in the region around \h\/\ J\\ ~ 2.00. As in the 
case with k = 0.60 the splitting of the isotropic-ordered 
transition into two transitions is only clearly observed for 
quite large system sizes L ~ 256, which prevents us from 
reaching conclusions about the nature of both transitions. 

B. Cumulant analysis of the order parameter 
distributions 

In the analysis of the phase transitions of model sys¬ 
tems it is quite useful to pay attention to the ratios be¬ 
tween the momenta of the order parameter distributions. 
Here we considered the ratios 34 =< O 4 > / < O 2 > 2 , 
i.e. 344 and 340 , which are closely related with the so- 
called Binder cumulants [46]. The analysis of the system 
size dependence of these quantities and, in particular, 
the crossings of the curves of 34 versus some thermody¬ 
namic field (T, h, ■ ■ ■) for different system sizes is often a 
very good choice to locate the phase transitions. In Fig¬ 
ure 9 we show the curves 34 (h) at constant temperatures 
for both order parameters and several system sizes. We 
can appreciate substantial qualitative differences in the 
shape of the 34 (h) functions defined on the order param¬ 
eter O t from the cases T/|Ji| = 0.80 (one transition), 
and T/|Ji| = 0.40 (two transitions with an intermediate 
nematic phase) and k = 0.6. For T/\Ji\ = 0.80 the tran¬ 
sition between the ordered phase (344 « 1 )to the disor¬ 
dered phase 344 —> 2 occurs quite abruptly as the system 
size increases, and the curves for different system sizes 
seem to cross at h/| Ji| ~ —0.305, which coincides (or it 


is quite close) to the crossing point of the corresponding 
lines for the 340 ratio. Notice that the maximum of the 
susceptibility \t (See Fig. 7) seems to happen exactly at 
the same value of h/| Ji|. In the case of T/\ Ji| = 0.40 the 
departure of 344 (h) from the ordered phase value, 344 = 1 , 
occurs at values of |h|/| Ji| clearly smaller than the corre¬ 
sponding departure of 340 (h). Looking at the cases with 
T/\Ji\ = 0.40, we can observe that the values of 340 at the 
crossings of the curves gi 0 (T) for different system sizes 
are consistent with the criticality of the two-dimensional 
Ising universality class [51], as one could expect from the 
symmetry of the order parameter O 0 . In addition, the 
results for the largest systems indicate that in the range 
of h/| Ji| values between the two transitions (as predicted 
by the maxima of the susceptibility \t) gAt(h) exhibits a 
plateau, with values consistent with 344 ~ 3 in the region 
where the nematic phase is supposed to be stable. The 
same type of results are found for the cases at k = 1 . 


C. Gallery of configurations 

In order to illustrate the differences between the stripe, 
nematic and disordered phases described in this work, in 
what follows we will present some representative configu¬ 
rations of the simulated systems for L = 256, considering 
the lattice gas version of the model. The following rules 
were applied to plot the configurations: (1) We consider 
only occupied sites (or cq = 1); (2) We plot segments 
between pairs of NN sites if, and only if, both sites are 
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h/IJjl h/IJjl 


FIG. 7. (Color online) Results for the J 1 -J 2 model with k, = 0.60 and T/|Ji| = 0.80. Different system sizes (as indicated in 
the legends) are considered. The results are shown as a function of the reduced external field, h/\Ji\. Panel a: Orientational 
order parameter (OOP); Panel b: Translational order parameter (TOP); Panel c: Magnetization per site, m; Panel d: \t = 
(dm/dh)T- 



FIG. 8. (Color online) Results for the J 1 -J 2 model with k = 0.60 and T/|Ji| = 0.6667. Different system sizes (as indicated in 
the legends) are considered. The results are shown as a function of the reduced external field, h/\Ji\. Panel a: Orientational 
order parameter (OOP); Panel b: Translational order parameter (TOP); Panel c: Magnetization per site, m; Panel d: \t = 
( dm/dh)r■ 


occupied; and (3) Four colors are considered, depend- direction depending on the value of the complementary 
ing on the direction of the bond (x or y), and for each coordinate. Each color is related with each of the four 
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h/IJjl h/IJ,l 

FIG. 9. (Color online) Order parameter cumulants, <74 as a function of the reduced external field for two reduced temperatures, 
different system sizes (see the labels and legends in the plots) and k — 0 . 6 . 


configurations in the ground state shown in Fig. 1. From 
this representation of the configurations, we expect for 
isotropic phases segments in both directions and four col¬ 
ors with similar probabilities; for nematic phases most of 
the segments will be in one of the directions and two col¬ 
ors will be predominant; whereas for the ordered phase 
most of the segments will have the same direction and the 
same color. We have chosen two cases. In the first one, 
shown in FIG. 10 we consider fixed value of the external 
field, h = 0 (k = 0 . 6 ) , and plot representative config¬ 
urations for three temperatures in the neighborhood of 
the transition temperature. It can be seen that there are 
no signatures of the presence of the Ising nematic phase. 
Above the transition temperature (left panel) one can ob¬ 
serve regions in the system where one of the four colors 
is predominant. At the estimated transition tempera¬ 
ture (middle panel in FIG. 10) the system has developed 
a large region with most of the bonds in vertical direc¬ 
tion and one predominant color, but still relatively large 
regions with the three remaining colors are still present. 
As the temperature is further reduced (right panel in 
FIG. 10) the regions with minoritary colors appear just 
as small islands. In the second case, shown in FIG. 11, we 
considered fixed temperature at T/\J\\ = 0.50 (k = 0.6). 
For this case we expect the stability of an intermedi¬ 
ate nematic phase. Configurations in the left (isotropic 
phase) and right (ordered phase) panels show similar fea¬ 
tures, apart from the lower density of segments, to those 
found in FIG. 10; however the configuration in the middle 
panel shows clear signatures of the nematic phase. Most 
of the bonds are oriented in vertical direction, but none 
of the two colors associated with this direction is predom¬ 
inant, and no long range order correlation in horizontal 
direction can be appreciated. 


VI. CONCLUSIONS 

We have shown that the J 1 -J 2 model in an exter¬ 
nal field has an intermediate phase with only orienta¬ 
tional order. This is the main result of the present work. 
We have performed an analytical approach based on the 
Cluster Variation Method in the square approximation 
and compared the results with Monte Carlo simulations. 
Both approaches are in qualitative agreement. We did 
not find evidence of orientational intermediate phases of 
nematic type for zero external field, where our results are 
compatible with those already known from the literature. 
Nevertheless, in the presence of an external field a phase 
with orientational but without positional order emerges. 
This is compatible with an Ising-nematic phase with Z 2 
symmetry, characterized in this context by the sponta¬ 
neous breaking of the Z 4 symmetry of the square lattice. 
We found that an Ising-nematic phase exists in a finite 
window of external field values and temperatures. 

For the parameter values studied we found that the 
disordered-nematic transition is of second order, with the 
order parameter going continuously to zero at ( h c ,T c ). 
Preliminary Monte Carlo results indicate that this tran¬ 
sition is probably in the Ising universality class. The 
nature of the second transition, from the nematic to a 
stripe phase with both orientational and positional or¬ 
ders, is more subtle. The CVM results give discontinuous 
transitions for the parameter values studied. Regarding 
the Monte Carlo results, it has been found that there 
are strong finite size effects. Then, simulations of very 
large system sizes are required to extract definitive con¬ 
clusions, which are beyond our present capabilities. The 
simulation results presented here strongly suggest the ex¬ 
istence of a stable Ising-nematic phase at low tempera¬ 
tures. In principle, according to the form of the order 
parameter O 0 , and to the crossings of the gi 0 {h) func- 
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FIG. 10. (Color online) Representative configurations for k = 0.6, h = 0, L = 256, close to the order disorder transition. From 
left to right T/|Ji| — 0.9756 (isotropic), T/|Ji| = 0.9718 (estimated transition temperature), and T/\Ji\ = 0.9662 (ordered 
phase with orientational and translational order). 



FIG. 11. (Color online) Representative configurations for k = 0.6, T/| Ji| = 0.50, L = 256, close to the order disorder transitions. 
From left to right h/\Ji\ = —0.400 (isotropic), h/\Ji\ = —0.397 (nematic phase), and h/\Ji\ = —0.390 (ordered phase with 
orientational and translational order). The typical length of the rods grows from left to the right. The densities in the lattice 
gas model are about (from left to right): 0.12; 0.21; and 0.42. 


tions for different system sizes, one expects a continu¬ 
ous transition from the disordered to the nematic phase 
with 2D Ising criticality. Regarding the transition be¬ 
tween nematic and fully ordered stripe phase it seems 
quite difficult to extract definitive conclusions with the 
type of calculations presented here. There are two ba¬ 
sic problems: first, at intermediate temperatures (say, 
T/|Ji| = 0.6667 for k = 0.6) the isotropic-nematic and 
nematic-stripe phase transitions are very close form each 
other, which makes difficult a finite-size scaling treat¬ 
ment. At lower temperatures, another difficulty arises. 
The correlation length (the length of the segments shown 
in the configuration plots) grows quickly as one reduces 
the temperature. This implies that, again, large systems 
have to be considered to extract conclusions. In addition, 
the nematic-stripe transition is clearly detected through 


the order parameter Ot , the functions ga or the max¬ 
ima of the susceptibility \t = ( dm/dk)T only for large 
system sizes. 
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